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ABSTRACT 

Using a hydrodynamic adaptive mesh refinement code, we simulate the growth and evolution of 
a galaxy, which could potentially host a supermassive black hole, within a cosmological volume. 
Reaching a dynamical range in excess of 10 million, the simulation follows the evolution of the gas 
structure from super-galactic scales all the way down to the outer edge of the accretion disk. Here, we 
focus on global instabilities in the self-gravitating, cold, turbulence-supported, molecular gas disk at 
the center of the model galaxy, which provide a natural mechanism for angular momentum transport 
down to sub-pc scales. The gas density profile follows a power-law cx r~ 8 / 3 , consistent with an analytic 
description of turbulence in a quasi-stationary circumnuclear disk. We analyze the properties of the 
disk which contribute to the instabilities, and investigate the significance of instability for the galaxy's 
evolution and the growth of a supermassive black hole at the center. 

Subject headings: galaxies: evolution — galaxies: high-redshift — galaxies: ISM — galaxies: nuclei — 
galaxies: structure 



1. INTRODUCTION AND MOTIVATION 

It is increasingly evident that the growth histories of 
different components of galaxies: stars, gas, and super- 
massive black holes, are intricately connected. Obser- 
vations indicate a relationship between the masses of 
supermassive black holes (SMBHs) and various prop- 
erties of their host galaxie s, such as the spheroid 
mass (M agorrian et al.l Il998f) and the velocity disper- 
sion of stars in the bulge fF crrarese fc Merrittl 120001: 
iGebhardt et al.l 120001 : iTremaine et all l2002ft . Recent 
simulations have shown that feedback from accreting 
SMBHs can regulate the growth of the black holes as well 
as the evolution of their host galaxies, making feedback 
a potentially important piece o f the theory of SMBH- 
host galaxy c o-evo lution (e.g. iDi Matteo et al.l 120051: 
Springel et al.ll2005t HJi Matteo et al.ll2007t ISiiacki et all 
2007ft. 

Building up the mass of a galaxy and driving AGN 
feedback requires a continuous replenishment of fuel in 
the center of the galaxy. A comprehensive understanding 
of such fueling is not possible without detailed knowl- 
edge about matter transport from large scales to the 
vicinity of the black hole. That transport helps deter- 
mine, in particular, both the amount of material avail- 
able for accretion, feedback, and star formation, as well 
as the strength and duration of the fueling. Large 
scale gravitational tidal fields created during major merg- 
ers of galaxies are thought to be effective at funneling 
matter toward the centers of galaxies. Indeed, simula- 
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tions and semi-analytic modeling of hierarchical growth 
scenarios have shown that a combination of accretion 
and black hole mergers c an effectively build-up the lo- 
cal black hole population dKauffmann fe Haehneltl | 200C 
Yoo fc Miralda-Escudel l2004t IVolonteri fc Reesl I2OOE 



Malbonetal.1 120061 : IVolonteri fc Reesl 120061: ILi et all 
2007ft . as well as assemble the ~ 10 9 M Q black holes 
already ob served to be prese nt in quasars at z « 6 
(|Fan et al.ll200l ILi et al.lL2007ft . 

Secular evolution of galaxies, particularly in the ab- 
sence of major merger events, can drive fuel down 
to the center as well. Global bar instabilities are 
thought to be efficient at transporting angular mo- 
mentum, allowing material in th e disk to move to- 
ward the center of the galaxy ([Roberts et all 19791 : 
c;™i,;« ~+ „i I h nonl. Im^„,,^u; h nod lci,i™„,.,„ „i I h nonl 



Simkin et al1ll980l: lNoTuchiill988t IShlosman et ailfl989l : 
Kormendv fc Kennicuttj l2004t iRegan fc Teubenl 12 004) . 
In the "bars within bars" scenario of IShlosman et al.l 
(1989), a large scale galactic bar drives gas inward where 
it forms a self-gravitating disk. As this disk becomes 
unstable, a smaller secondary bar forms, driving mate- 
rial down to scales of order 10 pc, at which point other 
physics can take over and transport the material the rest 
of th e way toward the black hole (e.g. IShlosma n et al.l 
[1990ft . 

Given the complexity of transporting matter from cos- 
mological scales all the way down to the circumnuclear 
region of a galaxy, and ultimately to the vicinity of a 
SMBH, it is clear that a thorough understanding of the 
relationship between SMBHs and galaxy formation re- 
quires modeling of a wide range of scales. Recently, 
cosmological SPH simulations have been combined with 
smaller scale simulations of a SMBH host galaxy environ- 
ment to study the effects of merger driven fueling and 
AGN feedbac k on supermassive black hole growth an d 
demography (|Di Matteo et al.ll2007t TSiiacki et al.ll2007ft . 
Reac hing a large dynamic range in N-body+SPH simula- 
tions, iMaver et all (|2007ft have studied galaxy dynamics 
during mergers resulting in a supermassive black hole 
binary. Small scale simulations have addressed gas dy- 
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namics in sub-galactic scale disks with high resolution 
(iFukuda et al.ll2000l IWac[^l20qil IWada fc Normanl l2001l: 
lEscalall2007t iWada fc Normanl 120071) .' finding the devel- 
opment of a turbulent, multi-phase interstellar medium 
(ISM). 

The study we present here follows the evolution of 
the circumnuclear region in a typical galaxy environ- 
ment rather than a dramatic, rare event such as a ma- 
jor merger, in a self-consistent cosmological simulation. 
We have specifically chosen a simulated galaxy that will 
evolve into a typical galaxy at z = 0. The simulation 
follows the galaxy during a phase of its evolution in which 
it has not undergone any major mergers within several 
dynamical time scales, in order to follow the development 
of instabilities in the circumnuclear disk, which may drive 
the transport of matter and angular momentum from 
large to small scales. The cosmological simulations use 
the adaptive mesh refinement (AMR) technique to self- 
consistently model the gas dynamics in a single galaxy 
at high resolution (sub-pc resolution in the center of the 
galaxy). A large dynamic range (> 10 million), achiev- 
able with AMR technique, allows us to bridge cosmologi- 
cal scales to scales relevant for molecular cloud formation 
(the birthplace of stars) and AGN fueling. After study- 
ing the physics in this basic model galaxy, we can begin 
to include physical processes that are directly relevant to 
the problem of SMBH growth in the context of galaxy 
evolution, such as AGN feedback. 

It is a complex task to implement mergers, feedback 
and secular evolution in large, cosmological simulations 
all at once. Our approach is to split the problem into 
pieces to be addressed one at a time, ultimately building 
a more realistic simulation. While the present paper fo- 
cuses on the structure of the galaxy, and the development 
of instabilities which may be responsible for driving mat- 
ter and angular momentum transport within the galaxy, 
we plan to explore the time evolution of such transport, 
specifically the evolution of the accretion rates of mass 
and angular momentum, in a subsequent paper. The na- 
ture of the circumnuclear region and a quantitative anal- 
ysis of accretion rates each have important consequences 
for galaxy evolution. 

The organization of the paper is as follows. In Sec- 
tion we describe the details of the simulation and the 
adopted "zoom-in" method. In Section^ we describe the 
features of the highly-resolved galaxy in a single zoom- 
in episode at z = 4, including a stability analysis of the 
disk and the potential role of physics not included in the 
current simulation. Section [4] examines the potential role 
of physics missing from the simulations. Finally, Section 
[5] contains a summary of our conclusions and their inter- 
pretation. 

2. SIMULATION 

Sections 12.11 and [2.21 briefly explain details of the simu- 
lation and of our methods. Additional details regarding 
tests of angular momentum conservation, and the treat- 
ment of dark matter particle discreteness effects will be 
available in the PhD thesis of the author Levine (also 
see the discussion of angular momentum conservation in 
Appendix [B]) . 

2.1. Simulation Code 



It is unfeasible to follow the evolution of a galaxy all 
the way down to the scale of a black hole hole accretion 
disk over cosmological times. Such a simulation spans a 
large range of scales (with a dynamic range > 10 7 ), and 
each spatial scale has a different temporal scale. As a 
result, the Courant condition for numerical stability re- 
quires very small time steps in the most highly refined 
regions, making it computationally expensive to follow 
the evolution of the galaxy over long periods of time. It 
is much more feasible to start with a lower resolution 
cosmological simulation, and zoom in to the small-scale 
region with increasingly smaller cell sizes and time steps. 
As the simulation evolves on small scales, the large scale 
portion of the simulation does not undergo much evolu- 
tion, and does not need as high resolution. 

The cosmological simulations presented in this pa- 
per are conduct ed with the Adaptive Refinement 
Tree CART) code dKravtsov et all [19971: lKravtsovlll999l ; 
iKravtsov et al.ll2002ft . The ART code includes the tech- 
nique of adaptive mesh refinement (AMR) , allowing high 
resolution of a galaxy residing in a small region of the 
cosmological simulation, while following the rest of the 
simulated region with lower resolution. The technique is 
appropriate for studying the structure and evolution of 
a galaxy over a large dynamical range and in a cosmo- 
logical context. 

The ART code includes a range of physics for modeling 
dark matter, stars, and gas dynamics. Gas is converted 
into stars in cells with densities greater than psf and 
temperatures less than T sp, where psf = 1.64M Q pc~ 3 
and T SF = 9000 K (see iKravtsovl 120031 for more de- 
tails), resulting in a star formation efficiency consistent 
with a Kennicutt law on kiloparsec scales (K ennicuttl 
119981) and with observations on 100 pc sc ales as well 
(e.g. lYoung eTaT]H99l iWong fc B lit zl [2001 . The code 
follows ISM physics, such as molecular hydrogen forma- 
tion, and gas cooling by heavy elements and dust under 
the assumption of collisional ionization equilibrium. The 
cooling and heating rates are tabulated as functions of 
gas density, temperature, metallicity, and redshift over 
the temperature ran ge 10 2 < T < 10 9 K using CLOUDY 
(jFerland et al.|[l998l ). which accounts for the metallicity 
of the gas, and formation of molecular hydrogen and cos- 
mic dust. In future studies, we plan to include the ART 
code's radiative transfer capabilities, which will be nec- 
essary for implementing AGN feedback. 

2.2. "Zooming-In" to the Center of a Galaxy 

We begin with a cosmological simulation, evolved from 
a realization of a random Gaussian density field at z = 
50, with periodic boundary conditions, measuring 6/i -1 
comoving Mpc across. The simulation was first run with 
low-resolution in order to select a galactic-mass halo for 
subsequent study. A Lagrangian region the size of five 
virial radii of the halo at z = was then identified at 
z = 50, and r e-sampled and run with higher resolution 
to z = 2.8 (see Kl ypin et al.ll200E for a description of the 
technique). The Lagrangian region of the simulation was 
automatically refined three levels, as a minimum. There 
are 2.64 x 10 6 dark matter particles in the Lagrangian 
region, each with mass 9.18 x 10 5 /i _1 M Q . Subsequent re- 
finement and de-refinement in this region followed a dark 
matter mass criterion, in which a cell was refined if its 
total dark matter mass is greater than 1.8 x 10 6 /i _1 M©. 



3 



At z = 4, the highest matter density peaks in the simula- 
tion have a maximum resolution of ss 50 pc (in physical 
units; corresponding to 9 levels of refinement on top of 
a 64 3 root grid). The largest halo has a total mass of 
2 x 10 11 M Q at z = 4, and it contains the progenitor of 
an L* spiral galaxy. This initial cosmological simulation 
was run including metal enrichment and energy feedback 
from supernova, and radiative transfer in addition to the 
physics described in Section |2~T1 

This particular z — 4 cosmological simulation was cho- 
sen for this study because it contains a galaxy that has 
not been disturbed by a major merger since it merged 
with a galaxy with 25% its mass at z rs 6. The time 
since this dynamically active period (~ 600 million years) 
is significantly longer than the dynamical time of the 
galactic disk at z « 4 (rs 15 million years). The galaxy 
is, however, far from quiescent, having grown in mass 
by w 25% from z — 5 to z = 4 via minor mergers and 
accretion from its environment. 

In the cosmological simulation, we allow further refine- 
ment (at z ~ 4), one level at a time, in a 1.5 kpc region 
centered on the galaxy, effectively zooming in to the cen- 
ter of the galaxy with increasing resolution. The "zoom- 
in" te chnique is similar to the one used by lAbel et all 
(2002) in simulations of the formation of the first star, 
which required an even larger dynamical range. Because 
we are starting simply in this first pilot study, and build- 
ing a more realistic simulation in pieces, radiative trans- 
fer and stellar feedback have been switched off in the 
zoom-in portion of the cosmological simulation (they will 
be examined in future studies). The high-resolution por- 
tion of the simulation employs additional refinement cri- 
teria, refining according to a level-dependent mass cri- 
terion on levels 11 and below 7 in the zoom-in region. 
The refinement criterion is defined so that the finer the 
resolution, the more aggressively the mesh refines, ensur- 
ing that there are enough highly refined cells to resolve 
structures on small scales. Specifically, the mesh refine- 
ment is super-Lagrangian in the circumnuclear region of 
the galaxy, and cells are marked for refinement if the gas 
mass in a cell is m lcvcl ~ 10 times the Lagrangian mass cri- 
terion (9 x 10 7 h~ M©), where m r = 0.7 is our fiducial 
value. The factor 0.7 was chosen through experimenta- 
tion, so as to make the refinement in the central region 
more aggressive, but not so aggressive that the simula- 
tion becomes too computationally expensive. 

During the zoom-in period, we increase the maximum 
level of refinement gradually, one level at a time, allow- 
ing the simulation to reach a quasi-stationary state on 
each level before moving to the next. The slow initial 
refinement allows the ART code to resolve spatial scales 
within the simulated galaxy while avoiding transient ef- 
fects. Since the time steps depend on the sound speed of 
the gas, they are therefore related to the dynamical time, 
and a fixed minimum number of time steps on each level 
forces the mesh to evolve for several dynamical times. 
We have run parallel simulations with 100, 300, and 1000 
minimum steps on each level, and have determined that 
a requirement of a minimum of 300 on each level suf- 

7 It should be noted that we adopt the convention of referring to 
level as the "top" level, and the maximum (most-refined) level as 
the "bottom" level, so that the terms "above" and "below" refer 
to lower and higher resolution, respectively. 



ficiently reduces transient numerical effects that result 
from refining too quickly. 

Throughout the initial zoom-in, we take precautions 
to avoid numerical artifacts. As the resolution of the 
simulation increases and throughout its subsequent evo- 
lution, the density of dark matter is smoothed on scales 
corresponding to level 9 resolution and above. This en- 
sures that each cell contains « 15 — 20 dark matter 
particles, avoiding effects resulting from their discrete- 
ness. It is well known that poor resolution in hydro- 
dynamics codes can lead to ar tificial fragmentation in 
the gas. iTruelove et al.1 ([1997) investigated these nu- 
merical instabilities and found that a simulation must 
resolve the Jeans length to avoid artificial fragmenta- 
tion. In addition to the mass criteria described above, 
the ART code's refinem e nt sch eme meets the Jeans con- 
dition of lTruelove et al.l (|1997[ ) on the maximum level of 
refinement (which increases gradually during the initial 
zoom-in episode), requiring that Ax/Aj < 0.25, where 
Ax is the resolution, or cell size, and Aj is the Jeans 
length. Additionally, as the mesh refines, rapidly collaps- 
ing structures can cause numerical instabilities at steep 
density gradients. The ART co de implements artificia l 
pressure support (as described in lMachacek et al] (120011 ) ) 
on the maximum level of refinement, in order to avoid the 
over-collapsing of structures. 

After the initial zoom-in, the highest level cells are 
0.03 pc across (corresponding to 20 levels of refine- 
ment), allowing us to determine the location for a su- 
permassive black hole test particle with high precision, 
but without resolving the black hole accretion disk, 
which would require additional physics (such as magneto- 
hydrodynamics, or MHD) not currently included in the 
ART code. Our results extend reliably down to a scale of 
at least 0.12 pc, or the size of four level 20 cells, which we 
consider to be our resolution limit. After reaching this 
maximum resolution (level 20), we replace 3 x 10 7 M© 
of the gas from the center of the simulated galaxy with 
a black hole point mass of equal mass and momentum. 
At present, we focus on the dynamical properties of the 
circumnuclear disk on scales where the gravity of the 
gas dominates that of the black hole. The addition of 
the black hole point mass here simply allows us to fol- 
low the black hole's location and velocity as a reference 
point. However the successful introduction of the black 
hole particle will be essential in subsequent simulations 
involving physics associated with the black hole (such as 
AGN feedback). 

After the initial zoom-in and the successful introduc- 
tion of the black hole test particle, the simulation then 
evolves at high resolution for approximately one dynam- 
ical time at the 100 pc scale, showing a highly resolved 
galactic disk. For a quasi-Keplerian disk, the number 
of orbital periods, iV or b, undergone by the simulation at 
radii less than 100 pc, is (i?/100pc)~ 3 / 2 . Figure Q] shows 
the 3-dimensional gas density in an approximately face- 
on view of the galactic disk, on several different spatial 
scales, for a single time step, demonstrating the large dy- 
namic range of the simulation. The distribution of stars 
is shown in Figure [H in both a face-on and an edge-on 
view of the galaxy on a scale of 3 kpc. The two dif- 
ferent colors correspond to old and young stars (yellow 
and blue, respectively), with the youngest stars in the 




Fig. 1. — Volume rendering of the 3-dimensional gas density on several different scales at z = 4, from 200 kpc in the upper left panel, to 
2 pc in the lower right panel. The color bar shows the density scale in units of cm -3 . Note the presence of the spiral arms and instabilities 
on a wide range of scales. (This figure is best viewed in color.) 



galaxy located primarily in the disk, and the oldest stars 
distributed more isotropically, populating the galaxy's 
bulge. 

3. RESULTS 

3.1. Structure of the Circumnuclear Disk 

Volume rendered images of the gas density in the z = 4 
simulation, such as those in Figure [TJ clearly show the 
spiral disk structure of the galaxy. The gas disk extends 
inward all the way to sub-pc scales, which are the small- 



est scales resolved in our simulation 8 . The geometric 
structure of the circumnuclear region is best illustrated 
by the eigenvalues of the inertia tensor, calculated in 
spherical shells around the black hole. The inertia ten- 
sor is given by 

8 In particular, the simulated disk does not contain the com- 
monly observed toroidal structure in the inner few pc of the galaxy. 
The failure of our simulation to reproduce the AGN torus may be 
the result of additional physics still missing in the simulation, as 
we will discuss in i| 1.11 
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Fig. 2. — The stellar component of the galaxy at z = 4 on a scale of 3 kpc. The stars in the image are divided into two populations, 
with yellow stars corresponding to the oldest stars in the galaxy, and blue stars corresponding to the youngest stars. The youngest stars 
are found primarily in the disk of the galaxy. (This figure is best viewed in color.) 
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Fig. 3. — Ratios of eigenvalues of the inertia tensor in spherical 

shells centered on the black hole. The ratio c/a <C b/a, indicating 0-1 I — — — — 

a disk structure. 0.1 1 10 100 10 3 
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I fr = m <xXjXk, (1) 

a=l 

where m a is the gas mass of cell a at radius r a = |x Q | 
inside a spherical shell containing N cells, and centered 
on the black hole. The eigenvalues of the inertia tensor 
define the principal axes of each shell, and their ratios 
describe the shape. Figure [3] shows the average radial 
profile of the ratios of the inertia tensor's eigenvalues. 
Above 0.1 pc, c is significantly smaller than b over several 
orders of magnitude in radius, indicating a disk structure. 
However, the fact that b/a is significantly less than one 
in parts of the disk indicates that the disk is not axially 
symmetric, and that it has large-scale structures, such as 
spiral waves and bars. 

The gas disk in the simulation is fully rotationally sup- 
ported, with the tangential component of velocity domi- 
nating over the radial component by at least an order of 
magnitude, as is illustrated in Figure 2J The top panel 
of Figure 2] shows the ratio of the radial to the tangential 
component of velocity, which is small throughout much of 
the disk, as the motion of the gas is almost entirely rota- 



FlG. 4. — Top: Ratio of radial and tangential velocity components 
of the gas, each averaged over 550,000 years. Bottom: Average tan- 
gential velocity component in units of the average quasi-Keplerian 
velocity. 

tional. The average radial velocity is negative, indicating 
the inflow of gas. The bottom panel shows the tangen- 
tial component of velocity in units of a quasi-Keplerian 
velocity, y/GM(r)/r, determined by the interior total 
mass, M(r), at each radius. The velocity is close to be- 
ing Keplerian, but since the mass distribution in the disk 
resembles a flattened ellipsoid, rather than a spherical 
distribution, the rotation is slightly super-Keplerian. 

A remarkable feature of the circumnuclear disk is that 
the average gas density profile, measured within spher- 
ical shells centered on the black hole, follows an almost 
perfect power-law with little evolution in time. Figure 
[5] shows the gas density profile from two different snap- 
shots of the simulation, as well as the average of the 
profile over a w 550,000 year period. Both the snap- 
shots and the average profiles of the gas density increase 
by w 8 orders of magnitude in the inner 100 pc of the 
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Fig. 5. — Radial profiles of the gas, dark matter, and stellar mass 
density. Snapshots of the gas density are shown at 100,000 and 
500,000 years after the initial refinement (thin solid lines; green and 
red respectively, in the color version), in addition to a 550,000 year 
average (thick solid line) . The snapshots do not show much vari- 
ation from the averaged profile, demonstrating that on timescales 
of several hundred thousand years, the disk is in a quasi-stationary 
state. The dotted line shows a power-law with slope —8/3 for com- 
parison, which matches the gas density slope well. The dark matter 
and stellar mass density profiles are shown by the long-dashed (blue 
in the color version) and short-dashed (magenta in the color ver- 
sion) lines, respectively. In the inner ~ 200 pc, the slope of the 
dark matter and stellar density profiles is ~ — 2 consistent with the 
adiabatic contraction model. 

simulation, obeying a steep power-law with slope —8/3. 
The stability of the gas density profile indicates that the 
disk is in a quasi-stationary state on timescales of sev- 
eral hundred thousand years. Figure [5] also shows the 
dark matter and stellar mass density profiles. In the 
central kiloparsec of the galaxy, gas comprises ~ 62% of 
the mass (~ 2.3 x lO 10 /^ 1 M total in the central kpc), 
dominating the dynamics of the disk. In contrast, the 
stellar and dark matter populations comprise ~ 13 and 
~ 25% of the galaxy's mass inside the central kilopar- 
sec. The disk is extremely gas rich at z = 4, because the 
galaxy is still actively growing and has not yet formed all 
of its stars. Interestingly, we find that both the stellar 
and dark matter profiles measured in the central ~ 200 
pc of the galaxy follow simple oc r~ 2 power-laws. The 
profiles match those predicted for the adi abatic contrac- 
tion o f dark matter from a n NFW profile (iNavarro et al.l 
119971) . using the model of iGnedin etaH (J2004|), all the 
way down to the 1 pc scale. 

Figure [5] shows that the RMS velocity dispersion of the 
gas (measured between neighboring cells, and therefore 
depending on the resolution) greatly exceeds the sound 
speed in the inner 100 pc, indicating supersonic turbu- 
lence in the disk. Supersonic turbulence decays on a dy- 
namical time scale, so the persistence of the turbulence 
in the circumnuclear disk of the galaxy indicates a driv- 
ing mechanism, which will be addressed in the following 
sections. The mean sound speed of the gas indicates a 
cold, molecular gas disk within the central 100 pc of the 
simulated galaxy. The floor in the sound speed shown 
in Figure [5] is determined by the minimum temperature 
to which our adopted cooling rates (computed using the 
CLOUDY package) apply. The mean sound speed is 
computed for cells with T < 20,000 K, because the spher- 
ical averages shown in Figure [6] include the low density 
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Fig. 6. — A comparison of the radial profiles of the tangential 
velocity component of the gas (solid), its RMS velocity dispersion, 
or turbulent velocity (dashed), and the mean sound speed, com- 
puted in cells with T < 20,000 K (dotted), time averaged over 
550,000 years. The sound speed is low, indicating cold, molecular 
gas within the central 20 pc. 



gas, infalling perpendicular to the plane of the galactic 
disk, and shock- heated to high temperatures (> 10 6 K). 
Additionally, turbulence inside the disk produces local- 
ized shocks that briefly heat individual cells to similarly 
high temperatures. The inclusion of shock-heated gas 
produces a broad temperature distribution and raises the 
mean temperature of the gas so that it does not effec- 
tively describe the typical sound speed of the gas in the 
disk. Therefore, the mean is computed over a tempera- 
ture range that selects cells in the disk, while excluding 
atypical cells, resulting in a more representative quantity 
for describing the thermal properties of the disk. 

3.2. Angular Momentum Transport 

The gas density profile of Figure [5] motivates a simple 
analytic description of angular momentum transport in 
the disk. Averaged over a sufficiently long time, the disk 
can be considered to be approximately azimuthally sym- 
metric. Additionally, the rotational velocity of the disk 
dominates over other velocity components: over the local 
velocity dispersion, and over the sound speed, as demon- 
strated in Figure [6] Therefore, the circumnuclear disk in 
the simulation is well approximated by a thin, rotation- 
ally supported, viscous, molecular gas disk. Conserva- 
tion of angular momentum in such a disk is desc ribed by 
the fo llowing equation in cylindrical coordinates (jPringld 
[Mil): 

S, T , 1 8 ,„ . 1 8G 
~di + HI)]?, = 27ri?^R' (2) 

where the angular momentum density normal to the disk 
is given by J z = _Ci? 2 _l, where _! = Vt/R. The surface 
density of the gas is S, and is the radial component 
of velocity, which is small compared to the rotational ve- 
locity of the gas. The right hand side of equation © 
describes the effects of a visco us torque, G, g enerated by 
turbulent motions of the gas. iPringld (|1981l ) parameter- 
izes the viscous torque G, which necessarily vanishes in 
solid body rotation, dVl/dR = 0, as 
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G{R 7 t)^2n^R 3 — (3) 

where v is the coefficient of turbulent viscosity. Dur- 
ing the initial zoom-in, the gas disk quickly reaches the 
power-law density profile shown in Figure [SJ On short 
time scales of a few hundred thousand years, traced by 
a single high-resolution zoom-in episode of our simula- 
tion, time-averaged radial motions through the disk are 
small compared to the rotational motion, as shown in 
Figure [H Over longer time scales, angular momentum is 
transported outward by a viscous torque driven by the 
turbulent motions of the gas, allowing accretion to slowly 
feed the black hole. The left hand side of equation J2]) 
can therefore be considered to be small on time scales of 
a few hundred thousand years. It then follows that the 
right hand side of equation should also be small, and 
that the torque G is approximately constant with radius. 

Given that the gas density p follows a power law (Fig- 
ure[5j), it is reasonable to approximate the surface density 
of the gas E by a power law with the radius R: 

E cx R~ p . (4) 

Since the gravitational potential $ satisfies V 2 <& = AirGp 
and E = 2pR, it follows that 

$ cx R 1 -?. (5) 

Since the disk is almost entirely rotationally supported, 
the tangential velocity x>t is effectively the circular ve- 
locity v c , defined by v 2 — R d^/dR. Thus, the angular- 
speed is 

and the viscous torque 

G cx R^ 1 -^ 2 !;. (7) 

If v cx R (which, as we will show in Section [3~3l is the case 
in our simulations), then equation ([7]) implies a power- 
law slope of (3 ~ 5/3. It then follows that the volume 
density of the gas should scale as p cx i? -8 / 3 , consistent 
with the slope measured in Figure[SJ Therefore, the sim- 
ple description of angular momentum transport, main- 
tained by turbulence in the disk, potentially provides a 
consistent interpretation of the power-law slope of the 
gas density profile of Figure GO 

3.3. Disk Stability and the Source of Turbulence 

In Section 13.11 we showed that the simulated galaxy 
contains a self-gravitating, turbulent, cold, molecular gas 
disk within the central w 100 pc. Such disks are suscep- 
tible to instabilities and to fragmentation. In the snap- 
shots shown in Figure the presence of instabilities is 
illustrated by waves moving through the simulated disk. 
The snapshots trace the evolution of disk structure on 
scales of w 125 pc, and w 12.5 pc (inset). The rota- 
tional period of the disk at the outer radius is such that 
the disk has undergone at least one full rotation (several 
more at smaller radii) over the time scales followed by 
the present simulation. The panels show the somewhat 
chaotic formation and re-formation of spiral structures 



on these scales, in contrast with the more ordered spiral 
structure seen on kiloparsec scales (Figure [1]) . This sec- 
tion includes an analysis of the behavior of structures in 
the disk on the scales shown in Figure [7] 

The Toomre Q parameter describes the stability of 
the disk against linear perturbations, and is given by 
Q = crturbK/7rGE, where k is the epicyclic frequency, 
and E i s the surface density of the gas, both determine d 
locally (|Toomrel [l96l iGoldreich &: Lvnden-BeH 119651) . 
Although the disk is cold, and the sound speed is low, 
the RMS velocity dispersion is substantially higher than 
the sound speed inside 100 pc, indicating that the gas is 
turbulent. Therefore, in place of the sound speed typi- 
cally used in the definition of the Toomre Q parameter, 
we have substituted a turbulent velocity, given by the 
quadrature sum of the sound speed and the RMS veloc- 
ity dispersion of the gas. The disk is, for the most part, 
quasi-Keplerian (as described in Section [3~Tj) . so that the 
epicyclic frequency k is proportional to the angular speed 
of the disk, f2 = v t /r. 

The top panel of Figure [5] shows the Toomre Q param- 
eter for the simulated disk, corresponding to the region 
shown in Figure [JJ Where Q < 1, the disk is susceptible 
to local gravitational instabilities resulting from axisym- 
metric perturbations. For Q > 1 the disk is likely to 
be stable against axisymmetric perturbations, but is still 
susceptible to instabilities arising from non-axisymmetric 
perturb ations, which are less sta ble than radial pertur- 
bations (jPolvachenko et alU QQ?). For the density profile 
shown in Figure [5J the disk becomes stable for Q > 2. 
The shaded region in Figure [H] (1 < Q < 2) therefore 
represents marginally stable values of the Q parameter, 
where the disk still might become unstable. Figure[5]sug- 
gests that the disk in the simulated galaxy lies mostly in 
the region of marginal stability inside ss 1 kpc. In the 
case of axisymmetric perturbations, the fastest growing 
unstable mode corresponds to the scale Af as t,r, at which 
dw 2 /dk = (where u> is the angular frequency of waves 
in the disk), given by 

2rr 2 

Afast.r - ~~Q^~- (°) 

Using the condi t ion for marginal stability from 
iPolvachenko et ail (|1997| ). the fastest growing mode for 
all modes (axisymmetric and non-axisymmetric), Af as t ia ii, 
is given by 

. Afast.r °"turb (n\ 
Afast.all — ^ ~~ ~GE~' 

The bottom panel of Figure [5] shows the ratio A/r for 
each of the fastest growing modes. In the inner 100 pc, 
the scales of the fastest growing modes are smaller than 
the radius by only a factor of a few, implying that the 
disk is stable on scales A <§C R (at least in the linear 
regime). This is an indication that perturbations in the 
disk operate on a range of scales, but always on scales 
that are an appreciable fraction of the size of the system, 
driving global instabilities, which generate turbulence on 
smaller scales. The disk remains locally stable all the way 
into a region less than 1 pc from the black hole. There is 
no catastrophic fragmentation into clumps small enough 
to form stars, so that angular momentum transport may 
continue uninterrupted by bursts of star formation. 
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Fig. 7. — The evolution of global instabilities in the disk over time, on a scale of m 125 pc and on a scale of fa 12.5 pc (inset). Each 
panel shows a volume rendering of the gas density (in cm -3 ) at a different epoch. The similarities between the appearance of the disk on 
the two different scales indicate the hierarchical structure of the disk. (This figure is best viewed in color.) 



By dimensional analysis, the turbulent kinematic vis- 
cosity v discussed in Section 13.21 can be described by 
Aoturb, where A is a characteristic scale for turbulence 
and ot ur b, ^ ne turbulent velocity shown in Figure [6l is 
a characteristic velocity. It is natural to identify the 
characteristic length scale with the length scale corre- 
sponding to the fastest growing unstable mode, Af as t,r, 
so that v ~ Afast^oturb- The top panel of Figure [9] shows 
that the turbulent viscosity v, given approximately by 
Afast.rOturb) scales linearly with the radius over much of 
the disk, so that v oc R. Although the behavior of the 
turbulent viscosity v can potentially change with time, 
the viscous torque G should remain approximately con- 
stant in radius for the above interpretation to be valid. 
The bottom panel of Figure [9] shows the average viscous 
torque G, normalized to its value at 10 pc. The normal- 
ized value is close to unity over most of the circumnuclcar 
region, which is indeed the condition needed to explain 
the —8/3 slope of the spherically-averaged gas density 
profile of Figure [5] 

Thus, a simple description of the structure of the cir- 
cumnuclear disk based on equations ((2][7|) and Fig. [9] pro- 
vides a good match to the simulation results, supporting 
the idea that global instabilities in the disk are responsi- 
ble for generating turbulence in the gas, resulting in the 
power-law slope described in Sections 13.11 and 13.21 
The conditions in the circumnuclear disk are poten- 



tially conducive to gaseous bar formation, the con- 
ditions for which have been studied extensively for 
different geometries and phy sical conditions, resulting 
in a variety of c riteria (e.g. lOstriker fc Peebles! 119731 : 
Efstathiou et alj _ ll982t IChristodoulou et all Il995al lbl; 
Bottemal 120031: IWvsell2004D. A com parison with the cri- 
terion of lOstriker fc Peebles! (|1973[ ), which characterizes 
stability in terms of the ratio t = T/\W\ of kinetic to 
gravitational potential energy, indicates that the disk is 
both secularly and dynamically unstable to bar forma- 
tion on all scales r > 0.1 pc adequately resolved by the 
simulation. Because of the chaotic and transient behav- 
ior of the instabilities, the structures shown in Figure [7] 
do not resemble a single, well defined bar, but rather a 
highly perturbed "bar." 

Since the Toomre criterion applies only to small lin- 
ear perturbations and the Ostriker-Peebles criterion de- 
scribes global modes, the question remains whether non- 
linear effects lead to fragmentation on smaller scales, be- 
low the resolution of the simulation. The top panel of 
Figure flOl shows the ratio of the local Jeans length of the 
gas, Aj, to cell size, Ax, for all simulation cells in the 
circumnuclear disk (with temperatures less than 10 3 K). 
In most cells, the Jeans length is larger than the cell size, 
preventing the gas from fragmenting into sub-cell sized 
clouds, ultimately leading to star formation. The Tru- 
elove criterion for preventing numerical fragmentation is 
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Fig. 8. — Top: Toomre Q parameter. For Q < 1 the disk is un- 
stable, and for Q > 2 the disk is stable. In the regime 1 < Q < 2, 
the disk is marginally stable (stable against axisymmetric modes, 
but not non-axisymmetric modes). Bottom: Average values of the 
the fastest growing unstable mode for axisymmetric perturbations, 
Af aB t r (solid), and the fastest growing unstable mode for all per- 
turbations, Af ast a n (dashed), each divided by radius. The ratio 
A/r > 0.1 throughout the circumnuclear disk, so there is no catas- 
trophic fragmentation down to small scales. 

enforced on the maximum level of refinement (level 20) 
only, so there are no numerical restrictions to prevent 
gas from collapsing all the way to level 20. However, 
only the central sub-pc part of the disk reaches levels 19 
and 20, as shown by the histogram in the lower panel 
of Figure [TOl indicating that the disk is stable to non- 
linear effects. While some cells in the disk have Aj < Air 
and may form stars, Figure [10] demonstrates that most 
of the disk is stable against collapse, and that there is 
no widespread burst of star formation. The few cells 
with A j / Ax < 1 may correspond to resonances where 
the pattern speed of waves in the disk is the same as the 
rotational velocity, leading to higher gas densities, and 
possible star formation. However, a detailed analysis of 
the behavior of these resonances with time is beyond the 
scope of this paper. 

3.4. The Nature of the Turbulence 

Turbulence has been widely studied in modeling of the 
ISMs of g alaxies and star forming regions (for recent re- 
views , see iMac Low fc K lcsscn 2004; iMcKee fc Ostriken 
2007). It is important to understand the nature of the 
turbulence in the present simulated galaxy as it plays 
a key role in the disk properties. The global instabili- 
ties arising in the simulated circumnuclear disk generate 
turbulence on a range of scales, maintaining the quasi- 
steady state of the disk. 

Turbulence dissipates energy in such a way that the 
turbulent velocity scales as uturb oc A 9 , where q = 1/3 for 



Fig. 9. — Top: Ratio of kinematic viscosity v, given approx- 
imately by Af ast r O" tur b, to radius R. The ratio v/R is roughly 
constant in radius. Bottom: The viscous torque G normalized to 
its 10 pc value. The viscous torque is constant in radius, consistent 
with the power-law gas density slope p <x R~ 8 / 3 . 
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Fig. 10. — Top: Scatter plot showing the ratio of the local Jeans 
length, Aj, to the cell size, Ax at 300,000 years after the initial 
refinement. Only cool cells with temperatures less than 10 3 K are 
shown. Bottom: Histogram showing the volume of level 20 and of 
level 19 cells as a function of radius. The simulation only refines 
to the maximum levels in the center of the circumnuclear disk. 

incompressible, sub-sonic (Kolmogorov) turbulence and 
5 = 1/2 for compressible turbulence in the zero-pressure 
limit (Burgers turbulence). The supersonic turbulence 
in the simulated disk falls between these two limits. In 
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Fig. 11. — Top: Turbulent velocity as a function of scale (or 
cell size) at 3 and 8 pc, for the fiducial and aggressive refinement 
runs. The more aggressive refinement run resolves the scaling of 
the turbulent down to smaller scales. Bottom: A comparison of 
the smallest resolved scale of turbulence and the fastest growing 
unstable mode, Af as t,all (where Af aat! aii is the same as in Figure[8j|. 



AMR simulations, it is straightforward to measure the 
turbulent velocity on different scales because the infor- 
mation on different levels of refinement is readily avail- 
able. However, the measurements are only accurate to 
within a factor of two in scale because the cell size de- 
creases by a factor of two with each refinement. The 
approximate scaling of the turbulence is shown in the 
top panel of Figure QT] The figure shows the turbu- 
lent velocity (given by the RMS velocity dispersion be- 
tween neighboring cells) at two different radii for two 
different simulation runs. The first is the fiducial run 
described in the previous sections, and the second is a 
short portion of the fiducial run, with a more aggres- 
sive refinement criterion on levels 13 and below, given by 
m ievei-io max jQ q.125] ( w here m r is an empiri- 

cally determined parameter defined in Section r2~2"]) . The 
more aggressive run probes smaller scales at a given ra- 
dius, in order to capture the scale of the turbulence. The 
turbulent velocities shown in Figure ITTI demonstrate the 
scaling of the turbulence down to small scales. It is dif- 
ficult to determine the precise slope of the scaling using 
measurements from the present simulation because of the 
limitation imposed by the refinement scheme. While a 
more precise characterization of the turbulence spectrum 
is beyond the scope of the present simulation, Figure 
ITTI shows that the slope approximately falls between the 
Kolmogorov and Burgers turbulence limits, as expected. 

At scales below the resolution at a given location in 
the disk, the turbulent velocities in Figure fTTI level off at 
constant values because the turbulence scaling can only 
be measured down to the resolution limit. The more ag- 
gressive refinement run (triangles in Figure lllj) probes 
smaller scales for a given radius, and therefore levels 
off at smaller scales. The bottom panel of Figure [Til 




Fig. 12. — A volume-weighted PDF of the gas density in cells at 
a scale of ~ 30 pc, with temperatures < 10 3 K (in order to exclude 
the hot corona outside the disk), averaged over ~ 550,000 years. 
The PDF is well fit by a log-normal distribution over the range 
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shows the approximate scale at which the flattening of 
the spectrum occurs for several different radii in the cir- 
cumnuclear disk of the fiducial run. For comparison, the 
fastest growing unstable mode, Af as t,aii is shown as well. 
The turbulent velocity measured at the resolution of the 
simulation is a numerical quantity, but it roughly corre- 
sponds to the turbulence on the scale of Af ast , r , which is a 
physical quantity. Figure [11] demonstrates that the sim- 
ulation resolves all scales that are Toomre Q unstable, 
supporting the resulting interpretation that instability 
driven turbulence maintains the quasi-stationary state 
of the circumnuclear disk. 

The density structure of the disk is also consistent 
with the description of turbulence in the disk. Super- 
sonic turbulence typically imposes a log-normal density 
distribution on the g demonst rated by models of 

turbu l ence in molecular clouds (e.g . Vazaucz-Sc madeiiil 
1994 iPassot fc Vazauez-Semadenil 119981 IWadal 120011 : 
Kritsuk et al.ll2007t iMcKee fc Os trikcr 2Q03), such as 
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where the mean, p, is defined as In p, and a is the dis- 
persion. Figure [12] shows the volume-weighted proba- 
bility distribution function (PDF), time-averaged over 
~ 250, 000 years, for a shell of gas at a radius of 30 pc. 
The PDF is well fit by a log-normal distribution over at 
least 4 orders of magnitude in density, consistent with 
models of supersonic turbulence. 

4. POSSIBLE EFFECTS OF MISSING PHYSICS 

The simulations presented here do not include all of 
the potentially important physics for galaxy evolution. 
Therefore, in this section we discuss the potential effect 
that the physics missing from the simulations might have 
on our results. 

4.1. Optically Thick Cooling 

Cooling in the high density, central region of the galaxy 
is not treated entirely correctly by the ART code. The 
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Fig. 13. — Column density of the gas. The horizontal lines show 
the column densities at which dust and H2 become optically thick 
to their own cooling radiation. 



column density of the molecular gas in the center is 
so high that the gas is expected t o become optically 
thick to its own cooling radiation (e.g. lRipamonti fc Abell 
2004). Additionally, the presence of dust grains in this 
region traps radiation and halts cooling furt her (see 
iDraine fc Leill984l : lOssenkopf fc Henningj[l994l ). Figure 
1131 shows the column density vertically through the cir- 
cumnuclcar gas disk. Horizontal lines are shown for com- 
parison, corresponding to the column densities at which 
dust and H2 each become optically thick to their own 
cooling radiation (~ 10 26 cm~ 2 and ~ 10 27 cm -2 , respec- 
tively, at 0.1 solar metallicity, which is the metallicity of 
the circumnuclear disk in the simulation at z = 4). In 
the inner rs 10 pc, the column density of the gas is large 
enough that the opacity of dust and molecular gas must 
be accounted for to accurately describe cooling in the 
simulation. 

The simulation runs presented here do not include op- 
tically thick radiative transfer. This may explain why 
the gas remains in a thin disk all the way in to sub- 
pc scales, and does not resemble the obscuring tori ob- 
served in the inner few parsecs of AGN. Should the opac- 
ity of the gas to its own cooling radiation be included, 
the gas around the mid-plane of the disk would not be 
able to cool and would heat to temperatures correspond- 
ing to the energy dumped into the gas by the turbulence. 
These higher temperatures may be able to provide a sub- 
stantial vertical pressure support in the central few par- 
secs of the disk, resulting in a thicker, more toroidal- 
like structure 9 . However, it has been suggested that 
hydromagnetic disk winds, and not hydrostatic pressure 
support, may be responsible for sustaining the optically 
and geometrically thick obscuration re gion, or "obscur- 
ing torus," in the nuclei of galax ies (e.g. lKonigl fc Kartid 
[1991 lElitzur fc Shlosmanl [20061 ). In which case, the in- 
clusion of optically thick radiative transfer may not be 
sufficient to create a "torus" in our simulations. 

In future simulations we plan to remedy this limitation 
of the present simulation by accounting for the opacity 



Fig. 14. — Equipartition (solid) and dynamical (dashed) mag- 
netic fields needed to influence the gas dynamics in the disk. 



of the high density region in the center, which will allow 
us to consistently incorporate radiative feedback from a 
central source and to test the above hypothesis. 

4.2. Magnetic Fields 

The ART code does not include magnetic fields, 
whereas MHD is certainly impo r tant for accretion disk 
physics (see iBalbus fc Hawlevl 119981 and references 
therein). It is for this reason that the present study 
has been restricted to scales larger than w 0.1 pc, cor- 
responding to about 10 4 Schwarzschild radii for a black 
hole of mass 3 x 10 7 M . 

At scales larger than the accretion disk, the absence of 
magnetic fields in the ART code is probably not impor- 
tant. An estimate of the strength that an equipartition 
magnetic field would need in order to affect the disk on 
small scales is given by i? 2 q « 47r ( oer 2 urb . In order to af- 
fect the large-scale dynamics of the gas, the magnetic 
field would have to have an even larger strength of or- 
der -B 2 yn ~ 47r ( ov 2 . Figure [14] shows estimates of B eq 
and -Bdyn, in order to demonstrate how high the mag- 
netic field would need to be to significantly influence gas 
dynamics in the simulated galaxy. For most of the galac- 
tic disk, the above estimates for the magnetic field are 
far high er than the few uG fields observed in real galax- 
ies (e.g. IZweibel fc Heilesl[l997l : lBeckll2001h . Even in the 
sub-pc region, where water maser observa tions indicate 
stron ger fields of a few tens of mG (e.g. iModiaz et al.l 
l2005t IVlemmings et al1l2007f) . the magnetic field is still 
too low to affect the gas dynamics in the model galaxy. 
We therefore conclude that magnetic fields, which are not 
included in our simulations, will have a negligible effect 
on the dynamical state of this disk unless their strength 
greatly exceeds observational measurements. 



9 Notice, that the outer, optically thin layers of such a disk would 
continue cooling efficiently, covering the hot interior with a cold 
molecular "skin" . Such a skin may become Rayleigh- Taylor unsta- 
ble, leading to fragmentation of the torus into individual clouds. 



5. DISCUSSION AND CONCLUSIONS 

We have used a large-dynamic range cosmological sim- 
ulation to study gas dynamics in the circumnuclear disk 
of typical mass galaxy at z sa 4 (evolving into an 
galaxy at z = 0), resolving the distribution of matter 
from megaparsec scales all the way down to sub-pc scales 
(with 20 levels of refinement). The simulation reveals 
a cold, fully molecular, self-gravitating, and turbulent 
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rotationally-supported gas disk, which is globally unsta- 
ble but locally stable. 

The global instability, operating on a range of scales 
comparable to the size of the system, generates turbu- 
lence down to the smallest scale resolved in the simu- 
lation. On small scales, the turbulence supports most 
of the disk against gravitational fragmentation and col- 
lapse. The disk, therefore, remains locally stable and 
reaches a quasi-stationary state. In that state, global in- 
stability drives bar-like and spiral- wave- like structures on 
time scales of the order of 100,000 years at 100 pc scales 
(and on proportionally shorter time scales at smaller 
radii), but on a 500,000 year time-scale the structure of 
the disk remains quasi-steady. 

The disk develops a power-law gas density profile with 
a well defined slope of —5/3 in surface density and —8/3 
in spherically-averaged volume density. This slope is a 
natural consequence of the dynamical state of the disk: 
the turbulence, generated by the global instability, re- 
distributes the angular momentum in the disk on a dy- 
namical time scale, reducing the gradient of the viscous 
torque. The resulting (quasi-)steady state uniquely de- 
termines the gas density profile we find in the simulation 
and the distribution of the angular momentum with gas 
mass. 

The turbulence in the disk drives the local outward 
transport of angular momentum and the inward flow of 
gas toward the supermassive black hole on time scales of 
at least 10 million years, which are too long to follow in a 
single zoom-in episode of the simulation. Thus, we only 
capture a single snapshot in the cosmological life of the 
supermassive black hole. In follow-up work, we plan to 
consider several snapshots taken at different cosmologi- 
cal times in order to describe the system on longer time 
scales. 

The dynamical state of the disk that we find in our sim- 
ulation appears to be consistent with the results of previ- 
ous simulations of isolat e d circumnuc l ear disks in galax- 
ies (iFukuda et all 120001: IWadal 120011: IWada fc Normanl 
I2nmb lEscalall2?KT7h . A distinctive feature of our approach 
is that we follow the dynamics of the circumnuclear disk 
within cosmological simulations. While most of the vol- 
ume of the simulation evolves little on time scales rel- 
evant for the dynamics of the circumnuclear disk, cos- 
mological scales provide realistic boundary conditions 
for the dynamics on sub-kpc scales. Since we find that 
the circumnuclear disk rapidly reaches a quasi-stationary 
state, its evolution is entirely governed by the boundary 
conditions. For the same reason, the fact that the cosmo- 
logical simulation that was used for the initial conditions 
did not resolve the scale of the circumnuclear disk, does 
not compromise our results: the quasi-stationary state of 
the disk does not depend on the initial conditions, and 
so the lack of power on scales below about 100 pc in the 
initial cosmological simulation is not important. 
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mological simulation to model the coalescence of two su- 
permassive black holes. While the detailed treatment of 
gas physics and the scientific questions answered by the 
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profile, although they do not elaborate on the physical 
mechanisms for the formation of this profile. 

The local stability of the disk, supported by highly su- 
personic turbulence, implies that the disk is capable of 
continuously feeding a central black hole, uninterrupted 
by catastrophic bursts of star formation, which could 
have consumed the available fuel were the disk locally 
unstable. This interesting dynamical state of the disk 
provides a potential solution to the problem of how AGN 
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not form an AGN torus on the appropriate scales (if it 
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mentation of gas cooling in the code. The ART code, like 
all existing cosmological codes, assumes that the cooling 
radiation from cosmic gas escapes freely into the IGM. 
This assumption breaks down at the densities and tem- 
peratures reached by the simulation in the inner 10 pc. 
At this scale the disk becomes optically thick to its own 
cooling radiation from dust and molecular hydrogen. In 
that regime the disk may to heat up and acquire a sub- 
stantial amount of pressure support, which may result in 
a puffier, more toroidal-like configuration for the inner 
several parsecs of the disk. A consistent treatment of a 
putative AGN torus will require simulations that incor- 
porate optically-thick cooling. 
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APPENDIX 
MEASURING THE VISCOSITY 

Here we address the validity of our estimate of the turbulent kinematic velocity parameter v, which motivates the 
interpretation that turbulent transport of angular momentum maintains local stability and the quasi-steady state of 
the circumnuclear disk. Instead of a direct measurement of the viscosity, we have assumed that the viscosity depends 
on characteristic velocity and length scales, given by a tm h and Af as t,r, respectively. The accuracy of this estimate can 
be tested by measuring the transport of gas due to turbulent motions. 

The ART code does not follow individual parcels of gas across cell boundaries, but rather advected quantities 
describing the gas. Therefore, in order to follow the turbulent motions of the gas, we introduce a passive scalar into 
the simulation which has a value of unity inside a spherical region centered on the black hole particle, and a value of 
zero outside this region, effectively "painting" the gas. As the simulation evolves, turbulent motions of the gas cause 
the profile of the scalar quantity to deviate from its original spherical form. The evolution of the profile depends on 
the turbulent kinematic viscosity, v, and in the early stages of diffusion, while the deviations from the initial profile 
are small, the profile approximately satisfies the linear diffusion equation for initial conditions, 

P(r,t = 0) = \ 1 > i ! r<ro ' (Al) 
[0, if r > r , 

(where tq is the radius of the initial sphere) allowing a more direct calculation of the viscosity. The solution to the 
linear diffusion equation for the above initial conditions, is 
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Fig. 15. — Mass-weighted profiles of the passive scalar as a function of radius, a nd e volving over time from a 30 pc sphere. The symbols 
are measurements from the simulation, and the lines are the best fits to Equation I A2I 



Figure [TBI shows samples of the profile evolving from a 30 pc sphere for different times (q). The measured profiles 
agree well with the solution given in Equation lA2l Starting from several painted regions at 1, 3, 10, 30, and 100 pc radii, 
we followed the evolution of the paint-weighted density profiles and measured the turbulent viscosity, v, as a function 
of radius. Figure \W\ shows the best-fit measurements of vjr in relation to the estimate given by CTt ur bAf as t ir /r. The 
lines in Figure [TBI show individual snapshots of o-turbAf as t,r A corresponding to the timescale over which the diffusion 
was followed (whereas Figure [HI showed a time-averaged estimate). The best- fit measurements match the estimates for 
vjr well, thus lending justification to the estimate for v used in the analytic arguments of Section [3.31 



Anomalous numerical transport of angular momentum is sometimes presented as a concern for adaptive mesh 
refinement simulations using interpolation schemes for calculating velocities on the mesh. We have tested conservation 
of angular momentum in the code by conducting a separate set of simulations of the collapse of an isothermal sphere 
(the details are included in the PhD thesis of the author Levine). The test was conducted for a mesh which refined 
along with the collapsing sphere, and for a pre-refined mesh, in order to test the conservation of angular momentum 
both as the mesh refines and as gas moves through the mesh. The sphere collapses into a thin disk [h/r << 1), which 
conserves angular momentum over several rotation periods. 

In Figure [T7] we demonstrate the conservation of angular momentum in the present simulation. As the maximum 
resolution of the simulation increases, the angular momentum profile initially evolves, as the viscous torque, G, 
described in Section 13.21 redistributes angular momentum within the disk. But after reaching the maximum, level 
20 resolution (the point at which we introduce the black hole particle into the simulation and let it evolve), the 
profile remains rather steady with time, because the disk has reached a quasi-stationary state. Figure [17] shows the 
angular momentum as a function of enclosed gas mass at several different times during a single zoom-in episode of 
the simulation. The figure shows the evolution of the angular momentum distribution from 200, 000 years before 
the introduction of the black hole particle, when the simulation has refined to level 11, to 500,000 years after the 
introduction of the black hole particle, when the simulation is fully refined, and has evolved for several hundred 
thousand years in a quasi-steady state. 

On the scale of the circumnuclear disk, the simulation has made several thousand time steps between t = 100 kyrs 
and t = 500 kyrs, meaning that the simulation has undergone significant evolution on this spatial scale. Numerical 
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FlG. 16. — Ratio of kinematic viscosity v, to radius R. The lines show estimates for u/R given approximately by Af aa t,rCturb! for three 
different snapshots. The points are estimates from a second method measuring the advection of gas at five different radii in the simulation. 
The snapshots correspond to the profile fits at different times, so the 1 pc point should be compared to the solid line, the next point to the 
dashed line, and the last three to the dotted line. 

effects, accumulating with each time step, would cause deviations in the angular momentum profile. Over longer 
time scales, the slow inward transport of matter will alter the angular momentum profile. However, on the hundred 
thousand year time scale we expect the profile to remain rather steady, reflecting the quasi-stationary state of the disk, 
unless there is anomalous numerical transport of angular momentum. Figure [17] shows no systematic deviations in the 
angular momentum profile after the initial refinement and the corresponding re-distribution of the angular momentum, 
demonstrating that there is no sign of unphysical numerical angular momentum transport on the time scales of the 
simulation. On time scales much longer than those simulated in the present paper, however, these effects may become 
important. 
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Fig. 17. — Angular momentum as a function of enclosed gas mass for several different times. The times are measured from the introduction 
of the black hole. The vertical dashed lines correspond to enclosed masses on scales of as 10 and »s 1000 pc. At each of these scales, the 
simulation has made approximately 230,000 and 3,600 steps, respectively (corresponding to the resolution at each scale), between t = 100 
and t = 500 kyrs. 



